Optical creation of vibrational intrinsic localized modes in anharmonic lattices 

with realistic interatomic potentials 
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Using an efficient optimal control scheme to determine the exciting fields, we theoretically demon- 
strate the optical creation of vibrational intrinsic localized modes (ILMs) in anharmonic perfect 
lattices with realistic interatomic potentials. For systems with finite size, we show that ILMs can be 
excited directly by applying a sequence of femtosecond visible laser pulses at THz repetition rates. 
For periodic lattices, ILMs can be created indirectly via decay of an unstable extended lattice mode 
which is excited optically either by a sequence of pulses as described above or by a single picosecond 
far-infrared laser pulse with linearly chirped frequency. In light of recent advances in experimental 
laser pulse shaping capabilities, the approach is experimentally promising. 
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I. INTRODUCTION 

Over the past several years, theoretical studies of the 
dynamics of anharmonic periodic lattices have estab- 
lished the existence of intriguing vibrational excitations, 
characterized by well-localized displacement patterns.u 
These so-called intrinsic localized modes (ILMs) can ex- 
ist at any site in a perfect lattice, in contrast to local- 
ized impurity modes in harmonic defect crystals. While 
vibrational ILMs with atomic scale localization have 
been obtained for increasingly realistic lattice-dynamical 
models,013 their experimental study has been impeded by 
the lack of direct methods for their excitation and veri- 
fication. Recently, ILMs in a complex quasi-lD charge- 
density wave system were inferred from resonance Ra- 
man data, through the use of a coupled electron- vibration 
model restricted to a single repeat unito 

Previously we have demonstrated theoretically that 
driven ILMs can exist as a steady-state response to an 
applied spatially homogeneous monochromatic driving 
force.Q As a natural extension of that work, and to ad- 
dress a key experimental question, we here describe the- 
oretically an avenue for the transient optical creation of 
ILMs: we show how they can be produced in a ID model 
lattice with realistic potentials by means of laser pulses 
whose time dependence is designed by an efficient opti- 
mal control scheme. 

Owing to their high power densities, lasers are attrac- 
tive sources for exciting large-amplitude ILMs. Indeed, 
there are some notable examples of experiments with 
powerful lasers in the regime of anharmonic lattice vi- 
brations. For instance, phonon resonances measured in 
the ferroelectric LiNb03 by experiments using single vis- 
ible laser pulses with a duration of 60 fs and an energy 
per pulse of 5 fiJ were interpreted in terms of "overtones" 
of the very, anharmonic lowest-energy TO phonon of A\ 
symmetry!] Also, experiments on Ti2 03 using single visi- 



ble pulses with 10 nJ energy and 70 fs widths to excite the 
Ai g mode apparently resulted in atomic displacements of 
about 0.07 A - corresponding to 2% of the interatomic 
spacing - and revealed anharmonic behaviorE 

Of more direct relevance here is the._fact that experi- 
mental laser pulse shaping techniquesErLi provide consid- 
erable flexibility in the time dependence of the applied 
force. Indeed, the use of tailored fields for vibrational ex- 
citation has attracted much interest, mainly in the con- 
text of optical control of dissociation and reactions in 
molecular chemistryjL9 but also for the selective expLLa=. 
tion of optical phonons in time-domain spectroscopy rW 2 \ 
We will focus on two methods for the transient opti- 
cal creation of ILMs: (i) ippulsive stimulated Raman 
scattering (ISRS) excitationEj by a sequence of femtosec- 
ond pulses at THz repetition rates from a laser operat- 
ing at visible or near-visible frequencies, and (ii) infrared 
(IR) excitation by a picosecond far-IR laser pulse. For 
both mechanisms, the system's dynamical response can 
be described classically, provided the underlying laser fre- 
quency for the ISRS case is well off resonance with vibra- 
tional and electronic transitions. 

Since ILMs are complex, large-amplitude vibrational 
excitations, the external fields necessary for their cre- 
ation are likely to have a complicated time dependence 
and large magnitudes. It is therefore advantageous to 
determine the optimal external fields by a systematic 
scheme. In engineering, the analogous task of designing 
the time dependence of an external force to steer a dy- 
namical system towards a desired target .state is a funda- 
mental problem. Optimal control theoryo provides a so- 
lution with a rigorous mathematical foundation, based on 
the variational minimization of a positive objective func- 
tional. In the realm of atomic dynamics, this approach 
has been successful in the design of external electric fields 
for selective. bond excitation in models of smalLiiarmonic 
latticescjiiZI and small anharmonic molecules.ta We ap- 
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ply a similar scheme to the creation of ILMs. 

The following section discusses our theoretical frame- 
work by providing the necessary details of the optimal 
control scheme and describing the specific anharmonic 
model lattice we use. Our numerical results concerning 
the optical creation of ILMs are given in Sec. Ill, which is 
divided into two parts according to the two different op- 
tical methods we consider. In Sec. IV we address aspects 
of the experimental feasibility and discuss our results. 
Section |V| concludes the paper, and an appendix pro- 
vides additional qualitative insight. Some o£-,the results 
presented here were summarized in a letter 



II. THEORETICAL BACKGROUND 

Although some notable exceptions have appearedjli 
most studies of ILMs have considered ID or 2D model 
lattices. This restriction to simpler model systems has 
facilitated progress towards a theoretical understanding 
of basic ILM properties, without the additional numerical 
complications .-encountered with more realistic 3D mod- 
els of crystals.B In addition to the reasonable expectation 
that many of these basic properties will transfer to the 
3D case, there is evidence that ID models may apply 
directly to some types of motion in 3D crystals. For in- 
stance, in Ref. || a realistic 3D lattice-dynamical model 
was considered and ILMs were obtained with displace- 
ment patterns localized along the edge of the crystal, 
which can be regarded as a natural generalization of a 
ID lattice. Furthermore, it is well known that along 
some high symmetry directions in 3D crystals, the har- 
monic lattice dynamics map onto an effective ID model 
involving the collective motion of lattice planes. We have 
performed preliminary studies showing that this mapping 
also occurs in the anharmonic case, for certain polariza- 
tion directions.E3 This point is addressed in more detail 
later. Hence, for a demonstration of the optical creation 
of ILMs, we consider a diatomic ID model lattice that 
incorporates realistic features, such as standard interpar- 
ticle potentials and measured harmonic properties of real 
crystals. 



A. The system: Hamiltonian and dynamics 

For longitudinal motion in an externally driven ID sys- 
tem, the Hamiltonian is 



H = J2 



pI 

2m„ 



V n,n-l 0„ - r„_ ; ) - f^\t)r n 



l>0 



(1) 



where particle n has mass m n , position r n and momen- 
tum p n and interacts with particle n — I via a potential 
Kn-i(r), to be specified below. The external force is 
given by ff*(t) = \V n £ 2 (t) and ff*(t) = q n £(t) for 



ISRS and IR excitation, respectively. Here £ (t) is the 
longitudinally polarized electric field, V n = (dV/dr n ) is 
the electronic polarizability derivative evaluated at the 
equilibrium configuration, and q n is the effective charge. 
With V n = {-l) n V and q n = {-l) n q, the external forces 
for both excitation methods have equal magnitudes and 
alternating signs: 



(2) 



The vibrational dynamics of this ID system are described 
by Hamilton's equations 



dH _ Pn_ 

dp n m n 



(3a) 



dH 



Pn = -^ r = /n({r m }) + (-l)"^) ) (3b) 



where 



fn({r m }) = [K,n-l( r n - r n -l)- 



l>0 



K+l,n( r n+l - r„)] 



(4) 



is the total internal force on particle n, with m (r n — 
r m ) = {dV n ^ m / dr)\ r=Tn ^ Tm . From Eqs. (^) we readily 
obtain the equations of motion 



InTn = fn({r m }) + (-l)".F(t), 



(5) 



which are more convenient for numerical simulations of 
the system dynamics. 



B. The scheme: optimal control theory 

In the following, we describe the most important as- 
pects of the optimal control scheme used in this work. 
More details can be found in Refs. [l5| and [ll| whose 
compact notation we adopt with some modifications. We 
consider a system of N particles. Unless explicitly stated 
otherwise, iV-dimensional vectors are denoted by lower 
case bold Roman letters and N x TV-dimensional matri- 
ces by upper case bold Roman letters. Moreover in phase 
space, 2iV-dimensional vectors are denoted by lower case 
bold Greek letters and 2iVx2iV-dimensional matrices by 
upper case bold Greek letters. 

We combine the TV-dimensional vectors r and p for the 
particle positions and momenta, respectively, to form a 
2TV-dimensional phase space vector 



(r 



T P T ) - 



(n, 



,r N ,Pi, ■ 



,Pn), 



(6) 



where the superscript T denotes the transpose. Equa- 
tions (0) arc then rewritten as 



£ = <f>[t,F(t)], 



(J) 
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with 



(8) 



where M is the Nx iV-dimcnsional diagonal mass matrix 
with elements (M)„; = m n S n i, (f)„ = /„({r m }), and 
the definition of the coupling vector q with components 
(q)„ = (—1)" allows a compact description of the exter- 
nal force terms. With this notation, we can rewrite Eq. 
(1 as 



M-r = f + qF(*). 
We furthermore specify initial conditions 
r(t = 0)=r h 

p(f = 0) = Pl 



(9) 



(10a) 



(10b) 



at t = 0. 

In order to apply optimal control methods to our prob- 
lem, we first need to define a positive functional that 
reflects the physical objectives to be reached. Starting 
from the lattice in some initial configuration (r-^pj), we 
want to excite a given anharmonic mode at a specified 
final time tf, while keeping the magnitude of the exter- 
nal force within reasonable limits. Combining the target 
mode positions and momenta pf into a final phase 
space vector ^, we define the objective functional 



^l>F j' dtf 2 (t), 



(11) 



where the nonzero elements (^) aa = ipa of the 2Nx2N- 
dimensional diagonal matrix ^ are positive weight fac- 
tors, as is ipp. Note that in order for all terms in the 
objective functional to have the same units, not all of 
the weight factors can be unitless. Furthermore, only 
trajectories £(t) satisfying the equations of motion (Q) 
are admissible during the optimization. Including this 
as a constraint in the objective functional, we obtain the 
modified objective functional 



J%F{t)] = J%T{t)] 



J tf dt\ T -{t-4>[t,r{t)]}, 

(12) 



where A T = (Ai, . . . , \2n) is a 2iV-dimensional vector 
of time-dependent Lagrange multipliers. This modified 
objective functional is minimized with respect to the tra- 
jectories £ and external force F(t) to obtain the optimal 
force. For clarity, these quantities are treated separately 
in the following two paragraphs. 

Variational minimization with respect to the trajec- 
tories {£a(i)}, including an integration by parts of the 



term f Q f dt A T • 5£ , yields dynamical equations for the 
Lagrange multipliers 



A + * ■ A = 0, 



(13) 



where $ is a 2iVx2./V-dimensional time-dependent ma- 
trix with elements {&) a p = 5(/>/3[£, .F(i)]/9£ a . At the 
final time tf, these dynamical equations are subject to 
boundary conditions 



\(t f ) = * • K(t/) - £/] 



(14) 



For our definition of JF(t)] [see Eq. (g)], the matrix 
$ decomposes into four ATxiV-dimensional blocks: 



* = 






K 


M- 1 






(15) 



where we have defined the iVx ^-dimensional time- 
dependent, symmetric dynamical matrix K with ele- 
ments (K) . = — df n ({r m })/dri. We can now simplify 
the description by considering separately the Lagrange 
multipliers (A r ) T = (Ai,...,Ajv) for the positions and 
(A P ) T = (Ajv+i, . . . , \2n) for the momenta. Then, Eq. 
( |l3| ) becomes 



A' - K • A p = 0, 



A p + MT 1 • A r = 0, 



which can be combined to yield 

M • X p = -K ■ A p . 



(16a) 



(16b) 



(17) 



We note that these dynamical equations for the momen- 
tum Lagrange multipliers correspond to the equations of 
motion for a lattice of N fictitious particles with masses 
m n and instantaneous positions A p , interacting via time- 
dependent harmonic forces. Combining Eqs. (|l4|) and 
dl6| ) , we obtain the corresponding boundary conditions 



A"(t / ) = ^.tp(t/)-p/] ) 



A*(ty) = -M- 1 ^ • [r(t f ) - r f ], 



(18a) 



(18b) 



where * p and \& r are A^xiV-dimensional diagonal weight 
factor matrices with elements (\& p ) n ; = 8 n iip n +M and 
{^ r )ni = 8 n iip n . Note, that these boundary conditions 
for the Lagrange multipliers are specified at the final time 
tf. Thus their determination requires the knowledge of 
the final positions r(tf) and momenta p(i/) of the par- 
ticles in the actual lattice. We emphasize that the intro- 
duction of Lagrange multipliers only serves the purpose 
of including the dynamics of the driven lattice as a con- 
straint in the objective functional (11). This results in 
coupled equations for the dynamics of the actual lattice 
and that of the Lagrange multipliers. 

The optimal external force is now obtained by minimiz- 
ing the modified objective functional (112) with respect 
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to T(t). This minimization can be done for an exter- 
nal forcep2 r (t) whose time dependence is allowed tejbe 
arbitrarjO or which has a prescribed analytic form.ta In 
view of the important aspect of experimental feasibility, 
we discuss the latter approach. For ISRS excitation, the 
laser frequency is neglected and we constrain T{t) to be 
a sequence of Gaussian pulses 



ISRS 



(*) 



(19) 



with individual heights {Si} and pulse center times {ti}, 
and common width A. For IR excitation, J-(t) is taken 
to have a linearly chirped IR frequency under a single 
Gaussian envelope with fixed width: 



T m (t) = Se {t -^ 2 /^ sin(0 + cut + at 2 ). 



(20) 



In both cases, the external force depends on a set of vari- 
able parameters {t^}, namely {Si},{ti} and A for the 
ISRS case, and S, to, 8, w, and a for the IR case. 

The gradient of the modified objective functional with 
respect to the external force parameters {r^} has compo- 
nents 



dt 



^(t,{ Tj }) + (X p ) T -q 



(21) 



where Eq. (||) has been used to write A T • d<p/dF(t) — 
(\ P ) T ■ q. The optimal control force F° pt (t) is obtained 
by finding the zero of this gradient. Following Ref. [ll], 
we use an iterative approach for this nontrivial numerical 
problem. Using an educated guess for the force param- 
eters {tj}, we integrate the dynamical equations (^) for 
the actual particles' driven motion forward in time from 
t = to t = tf, starting from the initial conditions (10). 
This yields K(i) for t in the interval [0, tf], as well 



as 



the final positions v(tf) and momenta p(tf), which are 
used to evaluate the boundary conditions ( |l8| ) for the 
Lagrange multipliers at t = tf. We next integrate Eq. 
( p7| ) backward in time from t — tf to t = 0. From this 
we obtain X p (t) for t in [0,i/], which is then used to 
evaluate the gradient components (|2l])- We adapted a 
fifth-order Gear predictor-corrector molecular dynamics 
(MD) methocO for the time evolution of the particles in 
the actual lattice and for the Lagrange multipliers. The 
force parameters are updated within a conjugate gradi- 
ent scheme,E3 and the procedure is repeated until the 
gradient ( pl| ) is zero to within a specified tolerance. 

For the case of purely harmonic potentials V n ^ n -i(r n — 
r n -i), 4> and K are time-independent matrices. Conse- 
quently, for initial conditions = and = corre- 
sponding to the lattice at rest in its equilibrium configu- 
ration, the exact optimal control force can be obtained by 
direct algebraic manipulation of the matrix equations (|7|) 
and (n6|) , without iteration. This is detailed in Refs. 



and [T^ for the application to harmonic molecules. Within 
this harmonic limit, we have used this algebraic approach 
as an independent test of our iterative MD method de- 
scribed above. 

The harmonic limit also allows considerations of the 
controllability of the system. A system is said to be con- 
trollable if an arbitrary specified target state £f can be 
reached exactly with a suitable external force T(t). In 
Ref. [It] it was pointed out that a harmonic system is com- 
pletely controllable only if all normal modes couple to the 
external force. Still, even if a system is not controllable 
in this rigorous sense, i.e., the external force couples only 
to a restricted set of normal modes, it may be possible 
to reach a final state that is close, although not 

exactly equal, to the specified target state £/. We will 
return to this aspect in Sec. Ill below. 



C. The model: interatomic potentials and 
characteristics 

We consider a ID diatomic lattice with masses to and 
M (> to), where nearest neighbors interact via Born- 
Mayer plus Coulomb (BMC) potentials 



V mM {r) = A mM e r/p 



-r/p _ <L 

) 

r 



(22a) 



(22b) 



while second neighbors interact via pure Coulomb poten- 
tials 



V mm (r) = V M m{t) 



'11 

r 



(23) 



and more distant neighbors are assumed to be non- 
interacting. Note that the interaction between an atom 
and its nearest neighbors distinguishes between "left" 
and "right" neighbors. Accordingly, minimization of the 
total potential energy of the static lattice leads to asym- 
metric nearest-neighbor equilibrium separations R r Q nM 
and R^ Im [^ i?™ M ). Although this asymmetry might ap- 
pear unusual at first sight, it indeed correctly represents 
the situation encountered by mapping the collective mo- 
tion of (111) planes in a zinc blende structure crystal 
onto an effective ID model lattice. E£l Furthermore, this 
asymmetry is in fact necessary to properly represent a 
lattice with first-order_B,aman active vibrational modes. 

In a detailed study,c3 we noted that ILMs should be 
classified according to their "associated" extended lat- 
tice mode (ExM), into which they spatially broaden with 
decreasing amplitude. Since we are focusing on the 
optical creation of ILMs, we study a model that ex- 
hibits ILMs associated with the optically- active ExM. 
This intuitively reasonable choice will later turn out 
to be essential. With an external force such as given 
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by Eq. (||), the anharmonic version of the optical 
zone-center mode (OZCM) with displacement pattern 
A{. . . , 1, —m/M, 1, —m/M, . . .) is both first-order Raman 
and IR active, where for the former the asymmetry 
R™ M ^ Rff m is necessary, as noted above. We there- 
fore need to find model parameters such that ILMs as- 
sociated with the OZCM exist. In Ref. ^3|, we obtained 
an ILM existence criterion based on the interplay be- 
tween fundamental harmonic and anharmonic dynamical 
properties of the ILM's associated ExM. For our case of 
interactions via realistic BMC potentials with their dom- 
inant soft anharmonicity, this criterion predicts the ex- 
istence of ILMs associated with the OZCM in diatomic 
lattices for which the optical branch of the harmonic dis- 
persion relation has a minimum at k = 0. In order to 
satisfy this constraint on the harmonic properties of our 
ID model lattice, while at the same time ensuring reason- 
ably realistic interatomic forces, we determined our BMC 
potential parameters by fitting the harmonic dispersion 
to branches of the measured dispersion for a real crys- 
tal along a high-symmetry direction. Since we consider 
both Raman and IR excitation, we focused on real di- 
atomic crystals with the zinc blende structure, for which 
the k = transverse optical (TO) phonon is both first- 
order Raman and IR active. In particular, the measured 
TO phonon branch of ZnS along the <111> direction 
exhibits the required minimum at k = (see Ref. 24). 
Using ZnS masses m — 32.1 amu and M = 65.4 amu, to- 
gether with A mM = 3.45 x 10 2 eV, X Mm = 2.73 x 10 3 eV, 
p = 0.279A, and q = 0.9e, we can approximately match 
those curves, as shown in the upper left panel of Fig. |l|. In 
the resulting model, R% lM = 1.67A and R$ Im = 2.95 A, 
and the maximum of the harmonic phonon gap occurs at 
the harmonic OZCM frequency ujq — 5.22 x 10~ 2 rad/ fs. 

re- 



= m n {2ufd n + - 



271 



cos(2^)/ n ({r m }), (25c) 



II B 



The optimal control scheme described in Sec. 
quires the positions and momenta of a specified target 
state. In order to obtain accurate predictions for the 
stationary solutions to Eq. (JsJ) without external driving, 
we use the well-established rotating wave approximation 
(RWA) for the particles' time dependence,!!] but general- 
ized to include static and second-harmonic contributions 
as well as oscillation at a mode's fundamental frequency 
to. The motion of particle n is assumed to be of the form 

r n (t) = b n + c n cos(ujt) + d n cos(2ujt) + r°. (24) 

After inserting this ansatz into Eq. (||), we multiply the 
resulting equations by either unity, cos(u>i), or cos(2wi), 
and average over a single period. For an iV-particle lat- 
tice, this yields a system of 3iV coupled nonlinear equa- 
tions for the static displacements {b n }, the fundamental 
dynamic displacements {c„}, and the second harmonic 
dynamic displacements {d n }: 



(25a) 



1 



2 k 



= m n u> 2 c n + - I d(f> coscp f n ({r m }), (25b) 
71" Jo 



where 4> — cot and / n ({r m }) is defined in (||). Once 
the boundary conditions are specified, these equations 
can be solved using standard numerical routines; the 
solutions in conjunction with Eq. ( |24| ) constitute the 
RWA. We verify our RWA predictions by performing 
direct MD simulations of Eq. (gL using a fifth-order 
Gear predictor-corrector methodEJ Traditionally, Born- 
von Karman periodic boundary conditions are employed 
for the description of bulk properties of "infinite" lattices. 
These boundary conditions are implemented by setting 
r Jl+ jv = fn + L, where L is the static-lattice equilib- 
rium length of the supercell. Following the nomenclature 
of Ref. |23|, we denote them standard periodic boundary 
condition (StdPBCs). While StdPBCs are convenient to 
describe infinite periodic lattices, we will also consider fi- 
nite systems with free-end boundary conditions (FBCs). 

Applying the RWA to our model lattice, we find that 
ILMs associated with the OZCM exist, as predicted by 
the criterion of Ref. ^3[ The ILM frequencies are in 
the harmonic phonon gap. The middle panel of Fig. 
HI gives the RWA static and fundamental dynamic dis- 
placements of such an OZCM-ILM at u> = 0.97o;o for 
a 22-particle system with FBCs. The displacements 
are normalized to the average equilibrium separation 
R = (R r nM + R$ Im )/2. For clarity, the corresponding 
second harmonic dynamic displacements, whose magni- 
tude is less than 8% of the fundamental dynamic dis- 
placements, are shown separately in the bottom panel. 
That this ILM is associated with a first-order Raman 
active ExM is reflected by the fact that its fundamental 
dynamic displacement pattern exhibits no inversion sym- 
metry. The static displacements seen at the ends of this 
finite system are the result of "surface" relaxation. The 
upper right panel shows the RWA frequency vs amplitude 
curves for the ILM and for the OZCM in the correspond- 
ing 40-particle lattice with StdPBCs. For each of these 
two modes the amplitude is given by the magnitude of 
the largest fundamental dynamic displacement c n . Also 
shown are the results of MD measurements of the ILM 
and OZCM frequencies, which agree to within less than 
0.5% with the RWA predictions. The to{A) curve for this 
ILM is indistinguishable from the corresponding curve 
for the 22-particlc FBC lattice of the middle panel, as 
one would expect from the mode's high localization. 

For our purpose of optical ILM excitation, the dynam- 
ical stability properties of the optically-active ExM are 
important, as will be seen later. We have therefore ex- 
amined the stability of the OZCM in our model lattice 
within an RWA-based approach detailed in Refs. ^ and 
p6| , but generalized to include perturbations of the sec- 
ond harmonic contributions. This stability analysis as- 
sumes infinitesimal displacement and velocity perturba- 
tions having an exponential time dependence exp(At), 
leading to a linear eigenvalu e pr oblem for the growth 
rates {A}. In previous studies,E3Ga we have shown that a 
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dynamical ExM instability with a purely real growth rate 
A is intimately connected with the existence of ILMs as- 
sociated with the ExM. The OZCM in our model lattice 
exhibits such an ILM-related instability. In the top panel 
of Fig. we show the predicted maximum growth rate as 
a function of the OZCM amplitude. MD measurements 
based on the "projection method" of Ref. |25| agree to 
within 8% with the RWA predictions, as indicated by the 
diamonds. By decomposing the instability perturbation 
into its spatial Fourier components, we can extract the 
wave vector (fc p ) max of the fastest-growing component, as 
discussed in Refs. ^3] and [25[ The bottom panel of Fig. || 
plots (fc p ) max as a function of the OZCM amplitude. At 
zero amplitude, (/c p ) max vanishes, and it increases to its 
maximum allowed value tt/(2Rq) over a restricted range 
of amplitudes. As discussed in Refs. ^ and|25|, the corre- 
sponding wavelength 27r/(fc p ) max introduces a preferred 
instability length scale at each amplitude. As the OZCM 
amplitude is increased from zero, the instability length 
scale decreases from infinity, reaching a maximum value 
of 4i?0: after which it remains constant for increasing am- 
plitude. Finite-time MD simulations of unstable OZCMs 
seeded with the fastest growing instability perturbation 
for various amplitudes reveal that the instability leads to 
a breakup of the OZCM into a periodic array of localized 
ILM-like excitations whose spacing is very close to the 
preferred instability length. T he sta bility properties of 
the OZCM will be used in Sec. |lIIB| below. 



III. OPTICAL CREATION OF ILMS 



Before we apply the optimal control scheme of Sec. 



[IB to the optical creation of OZCM- ILMs in the model 
lattice described in Sec. HC, we consider the question of 



controllability, which was briefly addressed in Sec. II B. 
We find that this aspect of our problem is different for 
different choices of the boundary conditions. 

For a harmonic diatomic lattice with StdPBCs, it is 
well known that the harmonic version of the OZCM is 
the only normal mode that couples to an optical-like 
force with a spatial dependence such as given in Eq. (|^). 
Hence, independent of the time-dependence of ^F{t), this 
external force can only excite the OZCM. On the other 
hand, if we start from zero initial conditions in an anhar- 
monic lattice with StdPBCs, this argument still applies 
at short times, since for the initially small amplitudes the 
interactions are dominated by the harmonic terms. It re- 
mains true even at longer times for our StdPBC model 
lattice, because the large-amplitude anharmonic OZCM 
continues to be an exact stationary solution. 

On the other hand, when a finite harmonic system with 
FBCs is considered, the external force (^) couples with 
appreciable strength to a set of normal modes. However, 
the size of this set decreases with increasing size of the 
system, and in the limit N — > oo of an infinite system it 
is again just the FBC normal mode corresponding to the 



OZCM in a StdPBC lattice which can be excited opti- 
cally. From these considerations of the harmonic case, we 
expect the external force (|J) to have better control over 
a system with FBCs than over a lattice with StdPBCs, 
when anharmonicity is included. However, with FBCs 
we also expect to see a dependence on the system size, 
with incre ased controllability for smaller systems. 

In Sec. Ill A we show how ILMs can be excited "di- 



rectly" in a finite system with FBCs. As expected, this 
approach depends on the size of the system. For the cre- 
ation of ILMs in crystal lattices, we exploit the fact that 
the OZCM is unstable and breaks up into ILM-like local- 
ized vibrations. Hence, although we can only excite the 
StdPBC OZCM directly, its decay can produce ILM-like 



vibrations "indirectly." This is detailed in Sec. Ill B 



A. "Direct" excitation of ILMs in finite systems 

1. ISRS excitation 

We first demonstrate the power of the optimal control 
scheme by showing our results for ISRS excitation of an 
OZCM-ILM in a 22-particle system with FBCs. Prelim- 
inary studies with simpler model systems revealed that 
the control scheme becomes more successful for longer 
control periods and larger numbers of pulses in the ap- 
plied sequence. Keeping the question of experimentally 
feasible laser pulse sequences in mind, we choose a control 
interval [0, tf = 50T ], where T = 2tt/ujo is the period of 
the harmonic OZCM. At t = the particles are at rest at 
their equilibrium positions. The system is then driven by 
a sequence of 49 Gaussian pulses. The MD time step used 
during the optimization for this case and for all of the 
simulations in this paper is To/100. For the target state 
at tf we specify the RWA-predicted positions and mo- 
menta of an ILM at frequency uo — 0.97wo- The displace- 
ment pattern is given in the lower two panels of Fig. [l]. 
Our initial studies with simpler systems also showed that 
a combination of weight factors tp n = l6(m n u>o/2) and 
4>n+N = 16/(m n o;o/2) for the target positions and mo- 
menta of particle n, respectively, and ip? = 1, guarantees 
a good balance of the various terms in the objective func- 
tional. Using these values in our control algorithmvields 
the pulse sequence given in the top panel of Fig. pi Th e 
common full width at half maximum (FWHM) 2vln2A 
of the pulses is 18 fs, and the amplitudes {Si} range 
from zero to 0.13 eV/A. The bottom panel shows that 
the application of this rather complex sequence of pulses 
in MD produces a strikingly "simple" result, namely the 
creation of a highly-localized excitation which persists al- 
most unchanged well after the applied field ends at t = tf. 

We note that the pulse sequence shown in the top panel 
of Fig. ^| does not represent a unique solution of the min- 
imization of the modified objective functional (|l2|). It 
turns out that depending on the initial guess for the pa- 
rameters of the sequence, the optimal control algorithm 
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reaches different local minima of ( |l2| ) with distinct val- 
ues. We have not studied this aspect in detail, but among 
the small sample of three different solutions we have ob- 
tained, the value of the minimum corresponding to the 
sequence shown in the top panel of Fig. ^ was the small- 
est. The application of each of the three sequences in 
MD simulations produced a long-lived stationary local- 
ized excitation like that shown in Fig. [|. From a practical 
viewpoint, the multiplicity of solutions suggests that it 
should be possible to impose additional constraints in 
the optimization algorithm to enhance the experimental 
feasibility of the resulting pulse sequence. 

Although the results of Fig. || appear promising, they 
are sensitive to the system's finite size as anticipated 
at the beginning of this section. For instance, apply- 
ing the identical pulse sequence to a 42-particle system 
with FBCs yields no localization. In Fig. ||, the system's 
ends are essential for the flow of vibrational energy to- 
wards the center to set up the target ILM. Hence this 
"direct" ILM excitation method is not suited for crys- 
tals, although it may have relevance for the excitation of 
anharmonic "local modes" in molecules such as benzene 
(CqHq). Indeed, ab initw-MD simulations for benzene 
readily yield local modes,EZl and the addition of an opti- 
mal control scheme may allow the prediction of optical 
waveforms for their creation. 



2. IR excitation 

We have attempted to obtain similar results for the di- 
rect IR excitation of ILMs in finite systems, using a single 
chirped far-IR pulse of fixed width for the external force. 
Although we can reach a final state that is quite 

well localized, the localization does not persist for any 
significant time after the applied field ends. Evidently, 
the more restricted prescribed analytic time-dependence 
of the chirped Gaussian far-IR pulse reduces the ability 
for successful control. 



B. "Indirect" excitation of ILMs in infinite lattices 



The results of Sec. [II A indicate that due to its size de- 



pendence, the direct excitation method is not a suitable 
approach for the creation of ILMs in crystals. However, 
our previous studies of the interrelation between ILMs 
and their associated ExM, as detailed in Refs. 23 and 
p5| , suggest an alternative approach: because we have 
designed our model such that gap ILMs are related to 
the optically-active OZCM, we can create localized vibra- 
tions in the gap, "indirectly" by optically driving the un- 
stable OZCM.E3 For this, our choice of a model in which 
the ILMs' associated ExM is optically active is essential. 



1. ISRS excitation 

To illustrate, we again use ISRS excitation. For the op- 
timal control algorithm target state, we specify an OZCM 
of frequency u> — 0.98w , m an 8-particle lattice with 
StdPBCs. This particular system size is large enough 
to avoid computational complications due to second- 
neighbor interactions across the supercell boundaries, 
and it is sufficiently small to accelerate the optimal con- 
trol scheme. The dynamics of the anharmonic OZCM are 
independent of the size of the StdPBC supercell, so that 
the resulting optimal fields apply to an infinite lattice. 
We choose this particular target O ZCM on the basis of 
our RWA stability analysis of Sec. II C. At ui — 0.98^0, 
the OZCM has an amplitude A = 4.37 x 10" 2 i? , 

a max- 
imum instability growth rate A max = 2.78 x 10~ 2 c<j, and 
the corresponding preferred perturbation wave vector is 
(fe p ) max = 0.6247r/(2i? )- Hence, this particular OZCM 
is expected to decay reasonably fast into ILM-like lo- 
calized excitations with an easily discernible spacing of 
about 6.4i?o- The advantage of choosing such an inter- 
mediate amplitude is clear from Fig. ||: at smaller am- 
plitudes the expected spacing is larger, but the corre- 
sponding growth rate is smaller, and vice versa for larger 
amplitudes. The issue of the size of the growth rate is 
important, because we want to ensure that the decay of 
the unstable OZCM occurs on time scales where other 
processes, e.g., damping, which are not included in our 
description of the lattice dynamics, will not significantly 
alter the results. However, since an amplitude corre- 
sponding to ~4% of the average equilibrium separation 
is quite large, we have repeated the optimization using 
as a target state an OZCM of smaller amplitude. The 
results for that case will be discussed at the end of this 
section. 

Starting from rest at t = 0, the system is driven with 
a sequence of 49 Gaussian pulses over a control interval 
of [0, 5 0Tq], as for the direct excitation detailed in Sec. 



Ill A 1| . We again use weight factors tp n — 16(m n u;o/2) 
and ip n +N = 16/(j7i n wo/2), but the simpler control task 
here allows us to increase the weight factor for the inte- 
grated square magnitude of the external force to ipp = 
10, without affecting the ability to reach the target state. 
In contrast to the pulse sequence given in Fig. |3J for the 
finite chain, the optimal sequence for OZCM excitation is 
found to consist of pulses having nearly equal amplitudes. 
Accordingly, we simplified the control algorithm so as to 
vary the pulses' common width, common amplitude, and 
individual pulse center times. 

The top panel of Fig. [| shows the resulting T{t) , which 
consists of pulses of FWHM 32 fs and amplitude 0.013 
eV/A. This is an order of magnitude less than the largest 
amplitude for the direct ILM excitation in the finite 22- 
particle lattice of Fig. |[ and the equal amplitudes render 
this sequence qualitatively simpler. However, a closer 
look reveals important details, demonstrating how the 
control algorithm has globally optimized this pulse se- 
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quence. The solid line in the upper panel of Fig. |3| plots 
the position of a light particle in the OZCM as a func- 
tion of time during the last fifth of the control interval. 
The thin vertical lines indicate the Gaussian pulse cen- 
ter times {ti} of the external force. As expected for an 
efficient impulsive driving force, the pulse centers {ti} co- 
incide with the zero crossings of the particle position. In 
addition, we measured the "instantaneous" amplitude of 
the OZCM during the excitation by taking half the dis- 
placement difference between adjacent turning points of 
the motion of a light particle and assigning it the time of 
the intermediate zero crossing. Then we calculated the 
corresponding undriven RWA frequency. This is shown 
as a solid line in the lower panel of Fig. ||, while the cir- 
cles indicate the frequencies 27r/(fct+i — ti) corresponding 
to the spacings between adjacent pulses of the external 
force. The good agreement between these two quantities 
brings out the important (and somewhat hidden) aspect 
of the sequence of Fig. |J: the spacing between adjacent 
pulses varies through the sequence in such a way as to 
maintain resonant impulsive driving of the anharmonic 
OZCM, whose frequency decreases as its amplitude in- 
creases over the control interval. In Appendix A we show 
that the qualitative aspects of this behavior can be read- 
ily understood by doing the optimization for a sequence 
of 5-function pulses. 

Applying our Gaussian pulse sequence to a 40-particlc 
StdPBC lattice in an MD simulation with the system 
initially at rest, we find that after the field ends at 
tf = 50T , the excited OZCM keeps vibrating with con- 
stant amplitude for about 150Tb until the perturbation 
due to accumulated computational round-off error trig- 
gers the instability and the OZCM decays into several 
localized excitations, as shown in the bottom panel of 
Fig. ||. The spatial array of localized excitations is not 
perfectly periodic because of the presence of instabil- 
ity perturbations of many wavevectors, with differing 
growth rates. Instead of relying on the clearly computer- 
dependent behavior of Fig. ^, we can provide the per- 
turbation necessary to trigger the OZCM instability by 
including the effects of nonzero temperature. The bot- 
tom panel of Fig. |^ displays an MD simulation for the 
same pulse sequence, but with random initial velocities 
corresponding to a lattice temperature of 5 K. The pres- 
ence of this perturbation triggers the OZCM decay much 
sooner. Of course the details of the MD results depend 
on the specific set of initial velocities, but for ten sets 
consistent with 5 K we find the same qualitative results 
as shown in Fig. |(| the ILM-like localized excitations 
resulting from the decay of the OZCM persist at fixed 
locations for several tens of vibrational periods and tend 
to move slowly through the lattice. 

As mentioned earlier, we have repeated the indirect 
ILM excitation via ISRS with the same control period, 
the same number of pulses in the sequence, and identical 
weight factors, but using as a target state a less anhar- 
monic OZCM at u> = 0.99^0, with the corresponding 
amplitude A = 3.03 x 10" 2 i? - In this case, our optimal 



control algorithm yields a sequence of pulses with 32 fs 
FWHM and amplitude 0.009 eV/A, compared with 32 fs 
and 0.013 eV/A for the target OZCM at u = 0.98w and 
A = 4.37 x 10~ 2 i?o. The qualitative time evolution of the 
unstable OZCM excited using this optimal sequence in 
MD simulations with nonzero initial temperature is simi- 
lar to that shown in Fig. ||, with the resulting ILM-like lo- 
calized vibrations having smaller amplitudes. Moreover, 
as expected from the larger preferred instability length 
scale and smaller growth rate at this OZCM amplitude 
(see Fig. |^), the localized excitations are further apart 
from each other and it takes roughly 20Xb longer than 
in Fig. |6| before a comparable degree of localization is 
reached. 



2. IR excitation 

We have also studied the indirect ILM creation via 
OZCM excitation using IR, assuming for T(t) a single 
Gaussian pulse [Eq. Q] of fixed width 33T (FWHM) 
and having a linearly chirped far-IR frequency over the 
control interval [0,100T ]- Using the same OZCM tar- 
get state and weight factors as detailed for the ISRS 
excitation above, the control scheme yields an optimal 
force with pulse amplitude 0.008 eV/A and chirp rate 
—2.0 x 10~ 7 fs~ 2 . The corresponding ^(t) is shown in 
the upper panel of Fig. [?], while the lower panel displays 
the results of applying this external force to a 40-particle 
StdPBC lattice in an MD simulation with initial random 
velocities corresponding to a lattice temperature of 5 K. 
Just as for the ISRS excitation of Fig. [| the external 
force excites a slightly perturbed OZCM, which subse- 
quently decays into ILM-like localized excitations. The 
above discussion concerning different sets of random ve- 
locities at the same temperature applies here as well. 

Again, the optimization was repeated with the same 
control period, the same fixed pulse width, and identical 
weight factors, but using as a target state the smaller- 
amplitude OZCM at uj = 0.99wo- In this case, our algo- 
rithm yields an optimal force with pulse amplitude 0.006 
eV/A and chirp rate —4.2 x 10~ 8 fs~ 2 . In analogy to 
the excitation of this OZCM with smaller amplitude via 
ISRS, the behavior in MD simulations with nonzero ini- 
tial temperature is qualitatively similar to that for the 
larger-amplitude OZCM shown in Fig. |7], but exhibits 
the same differences as in the case of ISRS excitation: the 
resulting localized vibrations are spatially further apart 
and it takes about 20To longer until a comparable degree 
of localization is reached. 



IV. DISCUSSION 
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A. Feasibility of the necessary external fields 

Having theoretically demonstrated the creation of 
ILMs using "designer" external forces, we now discuss 
the experimental feasibility of the corresponding fields. 
First we consider the excitation via ISRS. Pulse shap- 
ing for ultrashort (13 fs) visible laser pulses has been 
demonstrated^ with complicated final waveforms rang- 
ing from 12-pulse sequences with an overall Gaussian 
envelope and equal spacing to 6-pulse sequences with 
equal amplitudes and variable spacing. Furthermore, 
visible lasers producing ultrashort (18 fs) pulses with 
extremelji-jatge field magnitudes up to ~270 V/A are 
available. EJEII It is thus possible to produce the maxi- 
mum field strengths of 1.22 V/A and 0.38 V/A necessary 
for the examples of Figs. || and @, respectively, assuming 
a conservative value V — 2.5 A 2 for the polarizability 
derivative^ 2 ] However, the crucial experimental question 
is whether a given sample can tolerate such high electric 
fields in an experiment. .— ■ 

In a theoretical argument p3 Nelson and coworkers esti- 
mated the potential of single-pulse ISRS to excite large- 
amplitude anharmonic vibrations. Assuming pulses with 
10 /xJ energy focused to 50 fim. (FWHM) spot sizes, 
they predicted that a maximum phonon displacement of 
2 x 10~ 3 A could be produced in the organic molecu- 
lar crystal a-perylene. For pulse widths of 70 fs, these 
parameters correspond to a field strength of about 0.5 
V/A. Hence, their estimation of the displacement pro- 
duced by a single pulse and the field strength they consid- 
ered reasonable are comparable with the corresponding 
quantities in our indirect ILM excitation via ISRS. Fur- 
thermore, they suggested that in some materials coher- 
ent vibrational displacements in the 0.1-1 A range could 
possibly be achieved. These optimistic predictions were 
preceded by an actual experimental demonstration that 
single visible laser pulses of 70 fs width and 1 /iJ en- 
ergy focused to 150 /im spot sizes impulsively excited co- 
herent optic modes in a-perylene.O The field strengths, 
~ 0.05V/A, used in this single pulse experiment were 
one order of magnitude below those assumed for the sub- 
sequent theoretical prediction. However, as discussed 
in a later experimental paper by the same group|13 it 
turned out not to be feasible to use the proposed larger 
field strengths, because they exceed the fairly low laser- 
induced breakdown threshold of a-perylene. This se- 
quence of publications highlights how a material's laser- 
induced breakdown threshold can limit the potential of 
ISRS for the excitation of large-amplitude vibrations. 
However, it was pointed out in Ref. [l^ that materials hav- 
ing a substantially larger laser-induced damage threshold 
than a-perylene exist. 

The question of breakdown thresholds in ultrashort 
pulse laser-solid interaction is not yet very well studied. 
The breakdown thresholds for alkali halides under irradi- 
ation by near-visible laser pulses with pulse widths down 
to 10 ps were experimentally determined to be about 0.2 



V/A (Ref. |3^), with a tendency for the thresholds to 
increase with increasing frequency and decreasing pulse 
width. Only few measurements for femtosecond pulses, 
such as the pulses used in our ISRS excitation studies, 
are available in the literature. For fused silica and the al- 
kali fluorides, breakdown thresholds above 0.5 V/A were 
obtained with 275 fs and 400 fs pulses for visible and 
near-visible frequencies .Lj Values between 0.6 V/A and 
1.0 V/A were measured for fused silica, sapphire, magne- 
sium fluoride, and-glass, using visible laser pulses having 
a width of 120 fs.EZl Other experimental studies of fused 
silica obtained breakdown threshold fields of 3.8 V/A, 3.0 
V/A, and 3.3 V/A for visible-pulses of 150 fs, 100 fs, and 
55 fs width, respectivelyHI'Ea We did not find measure- 
ments of the breakdown threshold in ZnS in the relevant 
short-pulse regime. However, an experiment on optical 
coatings made from ZnS (Ref. ^) showed that its break- 
down field strength for nanosecond pulses is comparable 
to that of magnesium fluoride, for which the short-pulse 
values are given above. From these experimental results 
for pulses that are still 2-20 times wider than the ones we 



used for the excitation of ILMs via ISRS in Sees. Ill A and 



[II B, it appears that the maximum field strength for our 
direct ILM excitation of Fig. |3| for the 22-particle chain 
may be slightly too large to be realized in an experiment, 
while that for the indirect ILM excitation in the crystal 
lattice of Fig. M is below the values for breakdown. As 



discussed at the end of Sec. [HI B l| , the necessary external 
force magnitudes can be reduced by targeting an OZCM 
at a smaller amplitude, but this approach has limitations 
due to the competition between the time scales for the 
unstable OZCM decay and other processes in a crystal, 
e.g., damping. Moreover, the use of sequences with more 
pulses can also decrease the force magnitudes: accord- 
ing to the (5-pulse approximation (see Appendix A), the 
force magnitude is inversely proportional to the number 
of pulses in the sequence. However, since a larger num- 
ber of pulses requires longer control periods, the above 
caveat about the competition of time scales applies here 
as well. 

Turning to IR excitation, we note that the maximum 
force amplitude in Fig. [7] corresponds to a field strength 
0.008 V/A, assuming that q = l.Oe. Free-electron far-IR 
lasers produce picosecond pulses with intensities reported 
up to 4 x 10 7 W/cm 2 (Ref. |4l]), corresponding to a field 
magnitude of 0.002 V/A. Moreover, the frequency of free- 
electron far-IR lasers can be chirped at rates of —9 x 10 -9 
fs" 2 (Ref. g). These field magnitudes and chirp rates 
are within a factor of 5 and 20, respectively, of those 
used for the IR excitation of ILMs in Fig. 0. For this 
excitation mechanism, laser-induced breakdown should 
not play a role as a limiting factor, since the threshold of 
0.2 V/A obtained for alkali halides using near-IR pulses 
of 10 ps widthEj is well above the necessary maximum 
field m agnitud es obtained here. As discussed at the end 
of Sec. IIIB2, both the necessary field magnitudes and 
the chirp rates are reduced when a smaller amplitude 
OZCM is targeted. 
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We conclude that the fields necessary for the indirect 
excitation of ILMs via ISRS or far-IR as demonstrated 



in Sec. IIIB are reasonable and may be feasible in the 



near future. Among the approaches we have considered, 
excitation via ISRS seems more promising since lasers 
producing the necessary high field strengths are already 
available, although the experimental problem of laser- 
induced breakdown must be borne in mind. 

Preliminary results indicate that the indirect excita- 
tion of ILMs can also be achieved when simpler analytic 
time dependences for the external force are used, with- 
out necessitating significantly larger force magnitudes. 
For the indirect ILM excitation via ISRS, we repeated 
the optimization procedure using pulse sequences con- 
strained to have an equal variable spacing between the 
pulses as well as common variable widths and magnitudes 
and found that for the target OZCMs at u> = 0.98<x>o and 
lu = 0.99cjo the required force magnitudes were 2.3% and 
0.4% larger, respectively. Similarly, repeating the opti- 
mization for the indirect ILM excitation via IR using a 
single Gaussian pulse with unchirped variable frequency 
as well as fixed width and variable magnitude, we found 
the corresponding increases in the force magnitude to be 
5.0% and 1.3%. 

Another important experimental consideration is the 
robustness of the optimal fields. Considering the exper- 
imental limitations on the fidelity of shaped waveforms 
for ISRS, we note from Ref. ^9] that pre-specified pulse 
amplitudes and widths were reproduced to within 10% 
and pulse positions to within 10 fs. Randomly perturb- 
ing the ISRS pulse parameters of the sequence shown in 
the top panel of Fig. ^ for the indirect ILM excitation in 
the lattice within these margins reveals that although the 
perturbed excites the OZCM to a slightly different 
amplitude in each of the ten cases considered, a decay 
into localized excitations always occurs. On the other 
hand, the ISRS pulse sequence shown in the top panel 
of Fig. H for the direct ILM excitation in the 22-particle 
chain is more sensitive to such infidelities. In each of ten 
cases of randomly perturbing the ISRS pulse parameters 
within the above margins, a localized excitation at the fi- 
nal time results, but only in one case does this excitation 
persist after the applied field is turned off. Decreasing the 
margins to a 5% error for pulse amplitudes and widths, 
and 5 fs for the pulse positions, the success rate increases 
to four out of ten. 



B. Effects of nonzero initial temperature and 
damping 

We have obtained our optimal external forces assum- 
ing that the system is initially at rest, but in an experi- 
ment thermal fluctuations will always be present and we 
should consider their effect. For the case of indirect ILM 
excitation, we have seen that the efficacy of the external 
field is enhanced by the presence of velocity perturba- 
tions due to an initial temperature of 5 K, since they 



serve to trigger the OZCM instability. At this low tem- 
perature, thermal fluctuations are small enough to act 
just as perturbations on the zero temperature dynamics. 
Accordingly, application of the zero-temperature optimal 
force initially excites a slightly perturbed OZCM which 
subsequently decays, as seen in the lower panels of Figs. 
^ and 0, respectively. 

If we increase the initial temperature, the situation 
changes. We have performed MD simulations with the 
ISRS pulse sequence of Fig. ^ for ten different sets of ini- 
tial velocities appropriate to lattice temperatures of 77 K 
and 300 K. Already at 77 K the initial excitation can no 
longer be identified as a perturbed OZCM. Nevertheless, 
in all ten of the 77 K cases, the simulations result in well- 
localized ILM-like excitations, which appear sooner than 
for the 5 K case of Fig. |6[ However, fewer localized exci- 
tations are observed for the same size lattice, and an in- 
creased background of thermally excited long-wavelength 
acoustic vibrations is present. This trend continues as 
the initial temperature is raised to 300 K. Very similar- 
results are obtained when the optimal force for the far-IR 
excitation from Fig. |?j is used at higher initial tempera- 
tures. For comparison, we repeated the same nonzero ini- 
tial temperature simulations, but with no external force, 
and we observed no significant localization of vibrational 
energy. Therefore, although the optimal external force 
obtained for zero initial temperature does not achieve its 
original goal of exciting an unstable OZCM when used 
at these elevated temperatures, it nevertheless produces 
localized vibrations. 

Similar to the question of robustness with respect to in- 
fidelities in the pulse parameters, the optimal ISRS pulse 
sequence for the direct ILM excitation in the 22-particle 
chain is more sensitive to nonzero initial temperatures 
than is the sequence for indirect ILM excitation in the 
lattice. Performing MD simulations with the ISRS pulse 
sequence of Fig. || for ten different sets of initial velocities 
appropriate to lattice temperatures of 5 K, we find that 
a persisting localized excitation at the target site results 
in seven cases. Already at an initial temperature of 77 
K the success rate drops to zero, although in some cases 
a localized excitation is created at a site other than the 
target site. 

Vibrations in real crystals couple to other types of ex- 
citations and exhibit finite lifetimes - typically between 
5 and 300 phonon periods for optical phononsH3 We 
have included this aspect by adding phenomenological 
velocity-dependent damping to our MD simulations. An- 
harmonicity, which contributes significantly to phonon 
lifetimes, is already treated explicitly in our calculations; 
thus we assume a small damping constant corresponding 
to an OZCM lifetime of 100T . We repeated the opti- 
mization for the direct ILM excitation via ISRS of Fig. 
H after adding this damping, using the optimal pulse se- 
quence with zero damping as initial guess. This results in 
a qualitatively very similar sequence of pulses, but with 
a larger maximum force mag nitude of 0.24 eV/A, com- 
pared with the earlier result 0.13 eV/A for zero damp- 
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ing. Due to the presence of clamping the amplitude of 
the resulting localized excitation decays away within a 
few tens of To after the applied force ends. Similarly, 
when we include damping in the optimization of the in- 
direct ILM excitation via ISRS of Fig. [| we find that 
the common pulse amplitude increases from 0.013 eV/A 
to 0.016 eV/A. With T = OK initial conditions, the am- 
plitude of the resulting OZCM damps out before the ac- 
cumulated computational round-off error can trigger the 
OZCM instability, but for an initial temperature of 5 K, 
we again observe a break-up of the OZCM into ILM- 
likc localized excitations, whose amplitude subsequently 
damps out over few tens of OZCM periods. Additional 
calculations for the indirect ILM excitation via ISRS us- 
ing a damping constant appropriate to a shorter OZCM 
lifetime of 50Xb yield a common pulse magnitude 0.019 
eV/A along with qualitatively similar MD results. Thus 
we conclude that it is possible to create ILMs with damp- 
ing present, although the necessary force amplitudes are 
of course somewhat larger. Furthermore, detection of 
these ILMs would have to occur within few tens of OZCM 
periods after their creation. 

C. Relevance for real crystals 

While our ID model incorporates some realistic fea- 
tures, such as standard interparticle potentials and the 
measured harmonic dispersion of ZnS, it is not a model of 
any real crystal. We have also considered 3D models of 
ZnS-structure crystals using standard two-body central 
potentials between atoms out to second neighbors .Ej It is 
well known that the harmonic modes for k along <111> 
map onto an effective ID model involving the collective 
motion of (111) planes. We have shown that this map- 
ping also occurs in the anharmonic case, for certain po- 
larization directions. The resulting quasi-lD model, with 
effective anharmonic potentials between (111) planes un- 
dergoing collective motion with one transverse and one 
longitudinal degree of freedom, yields a representation 
of the actual motion in a 3D crystal and thus justifies 
our studies of purely ID latticescS Applying again the 
ILM existence criterion of Ref . we can choose the pa- 
rameters of the quasi-lD model such that the ILMs are 
associated with the optically-active ExM. Hence, we ad- 
justed our potential parameters so as to fit the measured 
harmonic phonon dispersion data of ZnS along <111>, 
as well as measured mode Griineisen parameters. Within 
the RWA this model exhibits gap ILMs associated with 
the anharmonic version of the k = TO phonon.Ej These 
results suggests that to the extent that the central poten- 
tials used in this model capture the anharmonic proper- 
ties of the real crystals, representative candidate materi- 
als for the indirect optical excitation of these ILMs would 
be ZnS, ZnSe, and the copper halides. We re-emphasize 
that these candidates for indirect optical creation of ILMs 
have two basic properties in common: (i) the k = TO 



phonon in these materials is both first-order Raman and 
IR active, and (ii) the frequency of the k = TO phonon 
is at the minimum of the TO branch along <111>, such 
that in conjunction with the fact that real potentials are 
dominated by soft anharmonicity, the criterion of Ref. |23| 
predicts the existence of gap ILMs associated with this 
k = TO phonon. 

V. CONCLUSION 

In conclusion, our optical excitation studies demon- 
strate theoretically that suitably tailored laser radiation 
offers a promising route for the laboratory creation of 
vibrational ILMs. The time dependence of the fields is 
determined by an efficient optimal control algorithm, de- 
signed to produce waveforms consistent with the rapidly 
developing experimental capabilities in laser pulse shap- 
ing. The direct excitation of ILMs in finite systems was 
demonstrated for impulsive stimulated Raman scattering 
by a sequence of ultrashort laser pulses at THz repetition 
rates, with variable spacing between consecutive pulses. 
For periodic lattices, ILM creation was achieved indi- 
rectly via decay of the unstable associated ExM which is 
excited optically either via multiple-pulse ISRS as above 
or via a single far-IR pulse with a linearly chirped fre- 
quency. For the indirect excitation approach, it is essen- 
tial to consider a lattice having ILMs that are associated 
with the optically-active ExM. 

The direct ILM excitation via ISRS requires rather 
complex pulse sequences and laser field strengths that 
may just exceed the breakdown threshold of a given sam- 
ple. Furthermore, this approach shows a dependence on 
the system size and is therefore not well suited for the 
creation of ILMs in crystal lattices, although it may be 
relevant for the excitation of local modes in molecules. 
Our studies show that the more advantageous means to 
excite ILMs in crystals is via the decay of their associated 
unstable anharmonic ExM, optically driven to a large 
amplitude. This approach not only requires smaller field 
strengths, but it succeeds even in the presence of thermal 
fluctuations and damping, and also with external forces 
constrained to have simpler analytic time dependences, 
i.e., pulse sequences with equal spacing between consec- 
utive pulses for ISRS and a single pulse with constant 
frequency for IR excitation. Hence, our studies point to 
a potentially fruitful avenue for experimentally accessing 
the regime of large-amplitude anharmonic vibrational dy- 
namics, which is very different than that for harmonic or 
weakly anharmonic systems. 
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J(K}) = E a ? 



X(E L -E f ), 



(A5) 



APPENDIX A: DELTA- APPROXIMATION FOR 
ISRS EXCITATION OF THE OZCM 

In this Appendix we demonstrate that the main char- 
acteristics of the optimal pulse sequence for the ISRS 
excitation of the OZCM in our model lattice described 



in Sec. 1IIB1 can be reproduced using a simplification 
which was considered in Ref. but in a different con- 
text. Here, the external force due to a sequence of short 
laser pulses is approximated by a sequence of L delta 
functions 



(Al) 



with individual positive amplitudes aj and pulse center 
times tj. In addition, we know that the dynamics of the 
OZCM in a diatomic lattice map onto that of an effec- 
tive anharmonic oscillator with mass m e e, displace men t 
x and momentum p. Including the external force (Al), 
the Hamiltonian for the driven case is given by 



tfcff 



p 



2m cff 



+ V e $(x)-xF(t), 



(A2) 



where V g(x) is an effective anharmonic potential. 

Starting with the effective oscillator initially at rest 
we want to excite it to a given final energy Ef while 
keeping the necessary external force magnitude minimal. 
This requires an optimization of t he parameters {aj} and 
{tj} for the external force ( |Al[ ). A single delta pulse 
a,j6(t—tj) boosts the momentum of the oscillator by aj at 
time tj and instantaneously changes the kinetic energy. 
Denoting the momentum of the oscillator immediately 
before and after the pulse by Pj-i and pj, respectively, 
we can write the energy boost as 



AEj 



2ajpj- 



2m 



2m 



off 



(A3) 



Evidently, for a given a,j the largest possible energy trans- 
fer occurs when the pulse arrives exactly at the time when 
Pj-i is positive and maximal, i.e., at the zero crossings of 
the displacement of the oscillator where its momentum 
is positive. Knowing the effective potential, this simple 
result is sufficient to determine the optimal pulse center 
times for a given set of pulse amplitudes {aj}. Further- 
more, we now know that the energy of the oscillator after 
the complete pulse sequence is given by 



Er 



We define an objective functional 




(A4) 



where the first term is just t he i ntegrated square mag- 
nitude of the external force (Al) and the second term 
constrains the energy El after the pulse sequence to be 
equal to the desired final energy Ef via introduction of a 
Lagrange parameter A. Minimizing this objective func- 
tional with respect to the pulse amplitudes {aj}, we find 
that the optimal pulse sequence consists of pulses with 
equal amplitudes 



1 



aj = —y / 2m c sEf. 



(A6) 



To summarize, the (^-approximation predicts that the 
optimal pulse sequence for the OZCM excitation con- 
sists of pulses with equal amplitudes whose pulse center 
times coincide with the zero crossings of the motion of the 
atoms. This agrees with our numerical results fr om the 



full optimal control scheme, as discussed in Sec. Ill B 1 
and illustrated in Fig. |5|. However, the ^-approximation 
underestimates the area under a single pulse by about 
20%, when compared with the pulses in the numerically 
obtained optimal sequence. Hence, this simplified de- 
scription correctly predicts the qualitative features of the 
optimal pulse sequence, but it cannot give quantitatively 
reliable results for the optimal pulse width and magni- 
tude. These require knowledge of the actual motion of 
the particles over the duration of each pulse and are ob- 
tained by applying the full optimal control scheme. 
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Roessler and Page, PRB, Figure 1 




particle 

FIG. 1. Characteristics of our ID diatomic BMC model lattice used for the optical creation of ILMs. Upper left panel: 
harmonic dispersion (solid line) and experimental transverse phonon frequencies (diamonds) along <111> in ZnS (Ref. 
Upper right panel: frequency vs amplitude curves for the optical zone-center mode (OZCM) (solid line) and the related intrinsic 
localized mode (OZCM-ILM) (dashed line) for 40 particles and standard periodic boundary conditions (StdPBCs). For the ILM 
curve the plotted amplitude is that of the mode's central particle, which is a light mass. The thin horizontal lines locate the top 
and bottom of the harmonic optical phonon band, indicated by the vertical bar. The circles and crosses are the results of MD 
measurements of mode frequencies for the OZCM and the OZCM-ILM, respectively. Our measurements and RWA predictions 
differ by 0.5% at most. Middle panel: static (squares) and fundamental dynamic (circles) displacements for an OZCM-ILM in 
a 22-particle diatomic BMC lattice with free-end boundary conditions (FBCs). The small (large) symbols represent the light 
(heavy) masses. The lower panel gives the corresponding second harmonic dynamic displacements as triangles. 
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Roessler and Page, PRB, Figure 2 



0.1 




FIG. 2. Stability properties of the OZCM in our ID diatomic BMC model lattice with StdPBCs. Upper panel: 
RWA-predicted maximum real instability growth rate of the OZCM as a function of the normalized amplitude (solid line). 
The diamonds give growth rate measurements obtained from MD simulations for a 40-particle lattice. Lower panel: wave 
vector of the fastest-growing Fourier component of the instability perturbation as a function of the normalized amplitude. To 
achieve good resolution, the stability analysis for both panels was based on a wave vector grid appropriate to a 2000-particle 
lattice. 
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Roessler and Page, PRB, Figure 3 
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FIG. 3. Direct ILM excitation via ISRS. Top panel: sequence of Gaussian pulses .F(t) for the direct creation of an ILM 
in a 22-particle system with free ends. The bottom panel shows the MD results of applying this sequence, with the particle 
displacements magnified by a factor 5. The applied field ends at tf. Note that the same force magnitude !F(i) acts on each 
particle. 
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Roessler and Page, PRB, Figure 4 
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FIG. 4. Indirect ILM excitation via ISRS with T=0 initial conditions. Top panel: sequence of Gaussian pulses T(t) for the 
indirect creation of ILMs in a 40-particle lattice with periodic boundary conditions. The bottom panel shows the resulting MD 
simulation, for zero initial conditions. Displacements are magnified by a factor 5, and only a portion of the lattice is shown. 
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FIG. 5. Details of the ISRS pulse sequence for the indirect ILM excitation. Top panel: position of a light particle during the 
OZCM excitation shown in the bottom panel of Fig. |^ as a function of time (solid line). Thin vertical lines denote the pulse 
center times t; of the corresponding optimal pulse sequence given in the top panel of Fig. [|. Bottom panel: RWA frequency 
calculated from the measured instantaneous OZCM amplitude during the excitation (solid line). Circles indicate the frequencies 
2ty/(U +1 -U). 
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Roessler and Page, PRB, Figure 6 



> 0.01 

CD 

CD 
O 

°0.00 



CO 

o 

Q_ 

Q) 
O 

c5 

Q_ 



0- 



— ' 



— .„ ~^W%ftwW«#J^^ 



+ 







T 



100 
time/T 



200 



o 



FIG. 6. Indirect ILM excitation via ISRS with T>0 initial conditions. Top panel: same as top panel of Fig. ^| The bottom 
panel shows the resulting MD simulation, for random initial velocities appropriate to a lattice temperature of 5 K. Displacements 
are magnified by a factor 5, and only a portion of the lattice is shown. 
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Roessler and Page, PRB, Figure 7 
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FIG. 7. Same as Fig. ^, but for the indirect ILM excitation by a single, linearly chirped far-IR pulse. 
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